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We study properties of excited states of an array of weakly coupled quasi-two-dimensional Bose 
condensates by using the hydrodynamic theory. The spectrum of the axial excited states strongly 
depends on the coupling among the various discrete radial modes in a given symmetry. By including 
mode-coupling within a given symmetry, the complete excitation spectrum of axial quasiparticles 
with various discrete radial nodes are presented. A single parameter which determines the strength 
of the mode coupling is identified. The excitation spectrum in the zero angular momentum sector 
can be observed by using the Bragg scattering experiments. 
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I. INTRODUCTION 

The experimental realization of optical lattices is 
stimulating new perspectives in the study of cold bosons. 
Optical lattices have enabled us to observe quantum phe- 
nomena such as number squeezing 0], collapses and re- 
vivals and the diffraction of matter waves (3| . Apart 
from these examples, BEG in optical lattices are particu- 
larly promising physical systems to study the superfluid 
properties of Bose gases |^ Q . The Bose- Hubbard model 
has been realized and the quantum phase transition from 
superfluid to a Mott insulator state was indeed observed 
experimentally 0, IE 13 ■ 1^ predicted that for deep 
optical lattices the condensate superflow can be lost not 
only by energetic instability but also by dynamical insta- 
bility 10, 11,1^]. The dynamic instability was verified by 
the experiments . In a seminal work by Kramer et al. 
|l4j . they have found the mass renormalization in pres- 
ence of the optical potential which decreases the value 
of the axial excitation frequencies. These discrete ax- 
ial excitation frequencies are experimentally verified 1151 . 
There are several theoretical calculations [1^, [l3 llM US 
for the sound velocity in a quasi-lD Bose gas placed in 
an ID optical lattice. 

The axial excitations of a cigar shaped condensate can 
be divided into two regimes: i) short wavelength excita- 
tions where wavelength is much smaller than the axial 
size, ii) long wavelength excitations where wavelength is 
equal or larger than the axial size of the system. In the 
former case, these excitations can be classified with a con- 
tinuous wave vector k. However, the finite transverse size 
of the condensate also produces a discreteness of the spec- 
trum. The short wavelength axial phonons with different 
number of discrete radial modes of a cigar shaped con- 
densates (without optical lattice) give rise to the multi- 
branch Bogoliubov spectrum (MBS) HIHI. The MBS 
was observed in a Bragg spectroscopy with a long dura- 
tion of the Bragg pulses >2^ . An array of weakly cou- 
pled quasi-two-dimensional condensates can be created 
by applying a relatively strong one-dimensional optical 
lattices to an ordinary three dimensional cigar shaped 



condensate along the symmetry axis. In presence of the 
periodic lattices, the MBS can be called as multibranch 
Bogoliubov-Bloch spectrum (MBBS). It is useful to study 
the MBBS in view of the possibility of the experimen- 
tal verification. It should be noted that all the modes 
in a given angular momentum sector are coupled among 
themselves for any finite value of the axial momentum. 
For example, when we excite the system to study the 
sound propagation along the symmetry axis, this per- 
turbation inherently excites all other low energy trans- 
verse modes having zero angular momentum. There- 
fore, all modes are coupled with each other in the same 
angular momentum sector. Martikainen and Stoof 
have studied the MBBS only for monopolc and lowest 
energy quadrupole modes without considering the cou- 
pling among the various modes within a given symmetry 
by means of time-dependent Gaussian variational ansatz. 
Later, Martikainen and Stoof 0| have calculated the 
spectrum of the phonon and the monopole modes by con- 
sidering only the coupling between the phonon and the 
breathing modes. But it is noted that the sound mode 
is coupled not only with the breathing mode but also 
with other low energy modes having zero angular mo- 
mentum. Similarly, the lowest energy quadrupole mode 
is also coupled with other low energy quadrupole modes. 
There is a lack of complete study on the MBBS in this 
system. For complete and correct description of MBBS 
we have to consider the couplings among all low energy 
modes in the same angular momentum sector. In our dis- 
cretize hydrodynamic description, the couplings among 
all the modes in the same angular momentum sector are 
included naturally and we will see in the next section. 

In this work we study the excitations in a stack of 
weakly coupled quasi-two-dimensional condensates. The 
multibranch Bogoliubov-Bloch spectrum of such system 
is presented by using the hydrodynamic theory. The 
MBBS strongly depends on the coupling between the in- 
homogeneous density in the radial plane and the den- 
sity modulation along the symmetry axis. Note that 
one can study only the spectrum of sound, monopole 
and quadrupole modes without considering the mode 
coupling completely by using the time-dependent Gaus- 
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sian variational method. Our discretize hydrodynamic 
method presented in this paper goes beyond the time- 
dependent variational method. In principle, we can cal- 
culate all low energy spectrum by including the mode 
coupling in a given angular momentum sector as long 
as the excitation energies are less than the chemical po- 
tential. We find that the multibranch Bogoliubov-Bloch 
spectrum changes due to presence of the mode-coupling 
within a given angular momentum symmetry. Therefore, 
the mode-coupling should be taken into account while 
calculating the spectrum correctly. 

This paper is organized as follows. In Sec. II, we con- 
sider an array of weakly coupled quasi-two-dimensional 
Bose condensates. Using the discretize hydrodynamic 
theory, we calculate the multibranch Bogoliubov-Bloch 
spectrum by including the mode coupling within a given 
symmetry. We give a brief summary and conclusions in 
Sec. III. 



II. MBBS OF A NON-ROTATING ARRAY OF 
BOSE CONDENSATES 

We assume that the bosonic atoms, at T = 0, are 
trapped by an external potential given by the sum of a 
harmonic trap and a stationary optical potential modu- 
lated along the z axis. The Gross-Pitaevskii energy func- 
tional can be written as 

= I dVi:Hr,z)[~^\/' + VUr,z) 

+ ||V(^,z)p + Kp(z)]V-(r,z). (1) 

Here, Vho{r,z) = ^{ujfr'^ + uj'^z'^) is the harmonic trap 
potential and Vop{z) — sE^ fiiv?{qz) is the optical poten- 

tial where E^ — is the recoil energy, s is the dimen- 
sionless parameter determining the laser intensity and q 
is the wave vector of the laser beam. Also, g = '^^^^ is 
the strength of the two-body interaction energy, where 
a is the two-body scattering length. We also assumed 
that ujr » so that it makes a long cigar shaped 
trap. The minima of the optical potential are located 
at the points Zj = jTr/q — jd, where d — ir/q is the 
lattice size along the z-axis. Around these minima, 
Vop{z) ^ M/2u1{z — ZjY , where the layer trap frequency 
is LOg = y/shq^ /M. In the usual experiments, the well 
trap frequency is larger than than the axial harmonic 
frequency, ujg » lo^. Therefore, we can also ignore the 
harmonic potential along the z-axis since the deep opti- 
cal lattice dominates over the harmonic potential along 
the z-axis. 

The strong laser intensity will give rise to an array 
of several quasi-two-dimensional condensates. Because 
of the quantum tunneling, the overlap between the wave 
functions between two consecutive layers can be sufficient 
to ensure full coherence. If the tunneling is too small, 
the strong phase fluctuations will destroy the coherence 



and lead to a new quantum state, namely Mott insulator 
state. 

In the presence of coherence among the layers it is 
natural to take the ansatz for the wave function as 

V'(x,y,z) = ^^j(a;,y)/(z - Zj). (2) 

Here, '(/'j (a^, v) is the wave function of the two-dimensional 
condensate at the site j and /(z — Zj) is a localized func- 
tion at j-th site. The localized function can be written 
as 

/(z-z,) = (M^)V4e-^(-^■)^ (3) 

Substituting the above ansatz into the energy functional 
and considering only the nearest-neighbor interactions, 
one can get the following energy functional: 

^0 = E / dxdv\-^^]-^i^,^v,,o\n^\ 
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+ ^E/ dxdyip^ip^ipjijj 

- J f dxdy['^j]^s4,,+^]^,+s]. (4) 

Here, J is the strength of the Josephson coupling between 
adjacent layers which is given as 

J = - j dzf{z)[-^Wl + Vop{z)]f{z + d) 

^ nu;r{^^r{^-4)se-^, (5) 

where = \J j^^^ ■ Also, the strength of the effec- 
tive on-site interaction energy is g^o — 9 J dz\fo{z)\'^ = 

f^(^), where Us = ^hjMujs. Eq. (g)) shows that 
each layer j is coupled with the nearest-neighbor layers 
J ± 1 through the tunneling energy J . The axial dimen- 
sion appears through the Josephson coupling between 
two adjacent layers. The Hamiltonian corresponding to 
the above energy functional is similar to an effective ID 
Bose-Hubbard Hamiltonian in which each lattice site is 
replaced by a layer with radial confinement. 

The Heisenberg equation of motion for the bosonic or- 
der parameter is 

(6) 

Using the phase-density representation of the bosonic 
field operator as 4'j = \f^j^^^^ ^^^^d neglecting the quan- 
tum pressure term, one can get the following equations 
of motion for the density and phase: 

- V'^jnj+isin(6'j+i -6'j)], (7) 
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and 



h0, - 



cos 



- Sj-i)] - Vho - .g2D? 



(8) 



Here, " • " represents the time derivative. In equilibrium, 
the condensate density at each layer is no(r) — 
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where we have neglected the effect of the tunneling energy 
J since it is very small in the deep optical lattice regime. 

Also, IjLq 



hujr^ y'S/iT^ is the chemical potential at 
each layer, where N is the number of atoms at each layer. 
In this system, we have two energy scales: the chemical 
potential of each layer fio s^^^ which is associated with 



the radial plane and the tunneling energy J 



se 



which is associated with the density modulation along 
the z-axis. The strength of the chemical potential can be 
enhanced by increasing the lattice depth or by increasing 
the number of atoms. The tunneling energy J decreases 
with the increasing of the lattice depth. 

We linearize the hydrodynamic equations around the 



equilibrium state, as rij = no 



and 9i 



The 



equations of motion for the density and phase fluctuation 
becomes 



2J 



—na{r)[2S9j - S9j-i - opj+ij 



d0i 



(9) 



and 



hS0j — —g2D^nj — 



J 



2no(r) 



[26nj — Srij^i — Srij+i]. (10) 



Note that the second term of the right hand side of Eq. 
(|10|l is proportional to the small parameter J and in- 
versely proportional to the large parameter no{r = 0) = 
/^o/ff2D- Therefore, we can neglect the term which is pro- 



portional to the 



J 



2no(r) 



After some algebra, we get second 



order equation of motion for the density fluctuation as 



- ^^^^no(r)[2(5nj- - (5nj_i - Jrij+i]. (11) 



The above equation tells us that the density fluctuation 
at each layer j is coupled with the nearest-neighbor layers 
j ± 1. We seek the normal mode solutions of the density 
fluctuations at layer j in the following form: 



5nj = 5n{r)e'''i^'^-'^'^^^^l 



(12) 



Here, k is Bloch wave vector (quasi-momentum) of the 
excitations. The Bloch wave vector p which is associated 



with the velocity of the condensate in the optical lattice 
is set to zero. 

Substituting the above equation into Eq. pi|l . we get 

—Lof{k)5n = ^^'^r-{no^rSn) ^^no{r) sm^{kd/2)Sn, 

(13) 

where I is a set of two quantum numbers: radial quantum 
number, and the angular quantum number, ni. The 
parameter J/i in front of the sin^(fcc?) term determines 
the strength of the coupling between the inhomogeneous 
density in the radial plane and the density modulation 
along the z axis. 

For k = 0. the solutions are known exactly and an- 
alytically [23. The energy spectrum and the normal- 
ized eigen functions, respectively, are given as, tuf = 
U!'^[\m\ + 2nr{nr + \m\ + 1)] and 

(14) 

Here, p!^'^'^ (x) is the Jacobi polynomial of order n and 
(f) is the polar angle. The radius of each condensate layer 
is Rq = 2iiq/Mlo1 and f = r/i?o is the dimensionless 
variable. 

The solution of Eq. (|13|) can be obtained for arbitrary 
value of k by numerical diagonalization. For fc 7^ 0, we 
can expand the density fluctuations as 



6n{r) = biSni{r, t 



(15) 



Substituting the above expansion into Ea. H13|) . we ob- 
tain, 

= [Cdf - [\m\ + 2nr{nr + |m| + 1)] 

- BQsin^{kd/2)]bi + Bo8in\kd/2)Y,Mivbi, 



where Coi — lui /cUr and the dimensionless parameter Bq is 
defined as 



Bo 



8J^o 



(17) 



The matrix element Mui is given by 

^^^^ ^ {l + 2nr+\m\) j J2~~2+\m\ + \m'\^^(m-m'),j, 

X p/^KI^°)(l-2f2)pa™l,o)(i_2f2). (18) 

The above eigenvalue problem is block diagonal with no 
overlap between the subspaces of different angular mo- 
mentum, so that the solutions to Ea. HlGI) can be ob- 
tained separately in each angular momentum subspace. 
We can obtain all low energy multibranch Bogoliubov- 
Bloch spectrum from Eq. (|16|l which is our main re- 
sult. Equations. (Qni and H18|l show that the spectrum 
depends on average over the radial coordinate and the 
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coupling among the modes within a given angular mo- 
mentum symmetry for any finite value of k . Particularly, 
the couplings among all other modes are important for 
large values of kd and Bq. It is interesting to note that 
the curvature of a mode spectrum depends on a single 
parameter Bq which is defined in Eq. (|17|l . The pa- 
rameter Bq can remain unchanged by changing values of 
the J and /zq in a various combination. Therefore, the 
curvatures of the spectrum of a given mode for various 
combinations of J and with fixed Bq are the same. 

Before presenting the exact numerical results, we make 
some approximation for a quantitative discussions. If we 
neglect the couplings among all other modes in the m = 
sector by setting I' — (^^,0) in Eqs. (fTB|) and l(TS)l. one 
can easily get following spectrum: 

cZ>2^ = 2nr{nr + 1) + (1 - A/„^,„jBo sm'^{kd/2). (19) 

The above equation can also be obtained by using first- 
order perturbation theory to Eq. Ijl^jl . In the limit of 
long wavelength, the rir — mode is phonon-like with 
a sound velocity Cq = V^^^fr, where AI* — is the 
effective mass of the atoms in the optical potential. This 
sound velocity exactly matches with the result obtained 
in Ref. and is similar to the result obtained without 
optical potential . This sound velocity is smaller by a 
factor of a/2 with respect to the sound velocity obtained 
previously for quasi-lD Bose gas placed 

in an optical potential. This is due to the effect of the 
average over the radial variable which can be seen from 
Eqs. lfTB|) and (P|l . 

In Fig.l we show few low-energy multibranch 
Bogoliubov-Bloch spectrum in the m = sector as a 
fmiction of kd by solving the matrix Eq. Hl()|) . 




kd 



FIG. 1: Plots of the low-energy Bogoliubov-Bloch modes in 
the m — sector. Here, J — O.lhuir and fio ~ 50hujr- Solid 
and Dashed lines are obtained from Eq. and Eq. p9^. 

respectively 

The lowest branch corresponds to the Bogoliubov- 
Bloch axial mode with no radial nodes. This mode has 
the usual form like ujr = Cgk at low momenta, where 



Cs is the real sound velocity. Note that Cg ^ cq which 
implies that the dispersion relations are modified due to 
the coupling among all other modes. The changes in the 
spectrum is clearly visible in the central part of the Bril- 
louin zone. This is due to the fact that the mode coupling 
is strong enough in the central part of the Brillouin zone 
due to the particular nature of the fc-dependcnt part (see 
Eq. (|16|l 'l. The second branch corresponds to one ra- 
dial node and starts at 2(jJr for fc = 0. The breathing 
mode has the free-particle dispersion relation and it can 
be written in terms of the effective mass (mj) of this 

mode as 1^2 (fc) — "^^r + ^7^- Fig.l shows that the mode- 
coupling does not affect on the breathing mode spectrum 
appreciably. The third and fourth lowest energy modes 
are also given in Fig. 1. These modes are also changed 
in the central part of the Brillouin zone due to the mode- 
coupling. One could sec from Fig.l that the effective 
masses of each modes are different. The group velocity 
along the z direction deviates from its long-wavelength 
limit when kd n. The mode coupling induced by the 
sin^(fco?) perturbation in Ea. l|ll|l becomes more signifi- 
cant with increasing fc and has the effect of lowering the 
sound speed. This coupling is associated with the inter- 
play of the density modulation along the z direction and 
the strong inhomogeneity of the equilibrium density in 
the radial direction in each plane. The effective masses 
are negative when kd > n. 

The coupling between the transverse quadrupole 
modes (m — ±2) and the modes in the m ^ ±2 sec- 
tor does not exist since these modes are orthogonal to 
each other as it can be seen from Eq. H18(l . However, the 
lowest energy quadrupole spectrum (n^ = 0,m = ±2) 
strongly depends on other low-energy quadrupole modes 
with various discrete radial nodes (n^ ^ 0,m = ±2). 
We neglect the couplings among all other modes in the 
m = ±2 sector by setting I' = {ur, 2) in Eqs. (fTH|l and 
(|18|l . then one can easily get following spectrum: 



Cbl =2 + 2nr{nr + 3) + (1 - M„,, 2-„,, 2)^0 sin^(fcrf/2). 

(20) 

In Fig. 2, we present first two low-energy MBBS for 
quadrupole modes. Fig. 2 clearly shows that the mode- 
coupling reduces the spectrum also for the quadrupole 
modes in the central part of the Brillouin zone. 

In Ref. j23|, the spectrum for the breathing and the 
lowest energy quadrupole modes are obtained analyti- 
cally within the Gaussian variational analysis. The mode 
coupling was not considered in this variational analysis 
[23]. In Fig. 3, we compare the spectrum of the breath- 
ing and lowest energy quadrupole modes obtained from 
Eq. H16() with those of obtained in Ref. . It is clear 
from Fig. 3 that the mode-coupling reduces the spec- 
trum strongly and it should be taken into account for 
calculating the spectrum correctly. 
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FIG. 2: Plots of the low-energy Bogoliubov-Bloch modes in 
the m = ±2 sector. Here, J — O.lhcJr and no = SOhur- Solid 
and Dashed lines are obtained from Eqs. HIGH and 1201 
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FIG. 3: Plots of the spectrum of breathing and lowest energy 
quadrupole modes. Here, J — O.lhuir and /lo = 507itJr. Solid 
and Dashed lines are obtained from Eq. 1161 and Ref. |2^. 
respectively 



III. SUMMARY AND CONCLUSIONS 

In this work, we have studied excitation energies of the 
axial quasiparticles with various discrete radial nodes of 
an array of weakly coupled quasi-two dimensional Bose 
condensates. Our discretize hydrodynamic description 



enables us to produce correctly all low-energy MBBS 
by including the mode couplings among different modes 
within the same angular momentum sector. We found 
that the mode-coupling strongly changes the spectrum. 
Therefore, it should be taken into account to calculate 
such spectrum correctly. The mode coupling is strong 
enough in the central part of the Brillouin zone. The 
single parameter Bq, defined in Eq. (|17|l . is identified 
which is always associated with the fc-dependent part 
and it scales with the product of two energy scales of 
this system, namely J and /io- The parameter Bq is a 
good measure for determination of the effect of the opti- 
cal lattices on the spectrum. Particularly, the spectrum 
for the phonon and breathing modes can be observed in 
a Bragg scattering experiments "26] as discussed below. 
The MBBS can be observed in the Bragg scattering ex- 
periments as the MBS was observed in Ref. ^22|. Due 
to the axial symmetry, the modes having only zero an- 
gular momentum can be excited in the Bragg scatter- 
ing experiments. In the Bragg spectroscopy, the conden- 
sate is excited by an external moving optical potential 
V = Vsit) cos{kz — cut), where VB{t) is the intensity of 
the Bragg pulses. This optical potential is created by 
using two Bragg pulses with approximately parallel po- 
larization, separated by an angle 9. The pulses have a 
frequency difference u; determined by two acousto-optic 
modulators. The wave-vector k is adjusted to be along 
the z-axis. Both the values of k and to can be tuned by 
changing the angle between two beams and varying their 
frequency difference. For small values of k the system is 
excited in the phonon regime and the response is detected 
by measuring the net momentum, Pz(a;, fc), imparted to 
the system by the Bragg pulses. The multibranch Bo- 
goliubov spectrum is obtained by observing the locations 
of the peaks in Pz{oJ, k) for various values of k. The fre- 
quency uj must be comparable to radial trap frequency 
LUr in order to excite the breathing and other modes. The 
duration of the Bragg pulses must be larger than the ra- 
dial trapping period, Tg > Itx jiOr in order to have large 
populations of the radial quasiparticle states. 



IV. ACKNOWLEDGMENTS 

This work is supported by a fellowship (P04311) of the 
Japan Society for the Promotion of Science. 



[1] B. P. Anderson and M. A. Kasevich, Science 282, 1686 
(1998). 

[2] C. Orzel, A. K. Tuchman, M. L. Fenselau, M. Yasuda, 
and M. A. Kasevich, Science 291, 2386 (2001). 

[3] M. Greiner, O. Mandel, T. W. Hansch, and I. Bloch, 
Nature (London) 419, 51 (2002). 

[4] Y. B. Ovchinnikov, J. H. MuUer, M. R. Dorey, E. J. D. 
Vredenbregt, K. Helmerson, S. R. Rolston, and W. D. 



Philips, Phys. Rev. Lett. 83, 284 (1999). 

[5] F. S. Catahotti, S. Burger, C. Fort, P. Maddaloni, F. 
Minardi, A. Trombettoni, A. Smerzi, and M. Inguscio, 
Science 293, 843 (2001). 

[6] S. Burger, F. S. Cataliotti, C. Fort, F. Minardi, M. In- 
guscio, M. L. Chiofalo, and M. P. Tosi, Phys. Rev. Lett. 
86, 4447 (2001). 

[7] D. Jaksch, C. Bruder, J. Cirac, C. W. Gardiner, and P. 



6 



ZoUor, Phys. Rev. Lett. 81, 3108 (1998). 
[8] D. van Oostcn, P. van dcr Stratcn, and H. T. C. Stoof, 

Phys. Rev. A 63, 053601 (2001). 
[9] M. Greiner, O. Mandel, T. Esslinger, T. W. Hansch, and 
I. Bloch, Nature (London) 415, 39 (2002). 
[10] B. Wu and Q. Niu, Phys. Rev. A 64, 061603 (2001). 
[11] V. V. Konotop and M. Salerno, Phys. Rev. A 65, 021602 
(2002). 

[12] A. Smcrzi, A. Trombcttoni, P. G. Kevredkidis, and A. R. 

Bishop, Phys. Rev. Lett. 89, 170402 (2002). 
[13] F. S. Catalioti, L. Fallani, F. Ferlaino, C. Fort, P. Mad- 

daloni, and M. Inguscio, New J. Phys. 5, 71 (2003); L. 

Fallani, L. De Sarlo, J. E. Lye, M. Modugno, R. Saers, 

C. Fort, and M. Inguscio, Phys. Rev. Lett. 93, 140406 

(2004). 

[14] M. Kramer, L. Pitaevskii, and S. Stringari, Phys. Rev. 

Lett. 88, 180404 (1998). 
[15] C. Fort, F. S. Cataliotti, L. Fallani, F. Ferlaino, P. Mad- 

daloni, and M. Inguscio, Phys. Rev. Lett. 90, 140405 

(1998). 

[16] K. Berg-Sorenson and K. Molmer, Phys. Rev. A 58, 1480 
(1998). 

[17] J. Javanainen, Phys. Rev. A 60, 4902 (1999). 



[18] M. Machholm, C. J. Pethick, and H. Smith, Phys. Rev. 

A 67, 053613 (2003). 
[19] E. Taylor and E. Zaremba, Phys. Rev. A 68, 053611 

(2003). 

[20] E. Zaremba, Phys. Rev. A 57, 518 (1998). 
[21] P. O. Fedichev and G. V. Shlyapnikov, Phys. Rev. A 63, 
045601 (2001). 

[22] J. Steinhauer, N. Katz, R. Orezi, N. Davidson, C. Tozzo, 
and F. Dalfovo, Phys. Rev. Lett. 90, 060404 (2003). 

[23] J. Martikainen and H. T. C. Stoof, Phys. Rev. A 68, 
013610 (2003). 

[24] J. Martikainen and H. T. C. Stoof, Phys. Rev. A 69, 
023608 (2004). 

[25] M. Fleesser, Andras Csordas, P. Szepfalusy, and R. Gra- 
ham, Phys. Rev. A 56, 2533(R) (1997); P. Ohberg, E. L. 
Surkov, 1. Tittonen, S. Stenholml, M. Wilkens, and G. 
V. Shlyapnikov, Phys. Rev. A 56, 3346(R) (1997). 

[26] J. Stenger, S. Inouye, A. P. Chikkatur, D. M. Stamper- 
Kurn, D. E. Pritchard, and W. Ketterle, Phys. Rev. Lett. 
82 , 4569 (1999); D. M. Stamper-Kurn, A. P. Chikkatur, 
A. Grlitz ,S. Inouye, S. Gupta, D. E. Pritchard, and W. 
Ketterle, Phys. Rev. Lett. 83 , 2876 (1999). 



